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1 .   Introduction 

The  purpose  of  this  paper  is  to  describe  two  methods 
for  constructing  smooth  bivariate  functions  which  take  on 
given  values  at  scattered  points  in  the  plane.   Given  the 
data  (x.f  y.,f.)/  i  -  1  ,...*!'»  both  of  these  schemes 
define  a  smooth  bivariate  interpolant   S  with  the  property 
that  S(x.,y.)  =  f.,  i  =  1,...,N. 

In  the  past  several  years  a  number  of  methods  have  been 
proposed  (e.g.,  [1],  [3],  [6],  [7],  [8],  and  [11]),  and  two 
survey  papers,  [2]  and  [10],  have  dealt  extensively  with  this 
topic.   The  problem  of  constructing  smooth  approximations 
based  upon  scattered  data  is  encountered  frequently  in  many 
areas  of  scientific  applications.   Common  examples  are: 
meteorological  information  such  as  rainfall  and  solar 
insolation,  geographical  information  such  as  elevations, 
geological  information  such  as  depths  of  underground 
formations,  and  engineering  data  such  as  stress  values 
obtained  by  finite  element  analysis.   A  somewhat  less  obvious 
example  is  given  in  [10]  where  it  is  described  how  human 
heart  potential  is  measured  at  irregularly  spaced  points  as 
an  aid  in  diagnosing  abnormal  heart  conditions. 

Since  a  number  of  methods  are  available  for  this 
important  problem,  a  project  was  undertaken  to  test  and 
compare  as  many  methods  as  possible  [4].   While  the  project 
included  both  global  methods  (meaning  that  S(x,y)  is  depend- 
ent on  all  data  points  regardless  of  their  distance  from 


(x,y))  and  local  methods  (S(x,y)  does  not  depend  on  data 
points  more  than  a  certain  distance  from  (x,y)),  for  large 
sets  of  data  it  is  necessary  to  use  local  methods.   During 
the  course  of  developing  and  testing  variations  of  previously 
published  schemes  we  have  discovered  two  which  appear  to  be 
preferable  as  general  purpose  methods.   While  neither  of 
these  methods  have  appeared  in  the  literature,  they  are  both 
modifications  of  previously  described  techniques.   In  Section 
2,  we  describe  the  basic  method  from  which  both  of  our  methods 
derive.   In  Section  3,  we  describe  the  details  of  the  two 
modifications  and  in  Section  4  we  show  the  results  of  applying 
these  methods  to  certain  test  problems.   In  Section  5,  we 
discuss  some  generalizations  and  ways  to  fine  tune  the 
schemes  to  suit  the  particular  needs  of  a  user. 

2 .   Inverse  Distance  Weighted  Least  Squares  Interpolation 
It  is  convenient  to  consider  the  data  as  coming  from 
an  underlying  function   f   and  to  view  the  interpolation 
process  as  an  operator  applied  to  this  function.   That  is, 
fi  =  f^xi'y±)'    i   =    lf-./N  and  S(x,y)  =  P[f](x,y).   Let 
p.(x,y)  denote  a  function  with  certain  properties  of  a 
distance  function.   In  particular,  we  assume  that  p.(x.,y.)  = 
0  and  that  1/p . (x.,y.)  is  a  nonnegative  decreasing  function 
as  (x,y)  "gets  further"  from  (x.,y.).   For  example,  the 
usual  Euclidean  distance,  d.(x,y)  =  V(x  -  x-)   +  (y  -  y-) 
is  one  possibility  for  p..   We  denote  by  cj)  . ,  j  =  l,...,m 
a  set  of  basis  functions  to  be  used  for  a  least  squares 
approximation . 


McLain  [7]  was  the  first  to  consider  a  family  of  inverse 
distance  weighted  least  squares  approximations.   The  general 
form  of  these  interpolants  is 


m 
(2.1)      P[f](x,y)  =  I      a.  (x,y)cj).  (x,y) 


j-l   ^ 


where  a.(x,y),  j  =  1, ...,m  represents  the  solution  of 


N 
(2.2)      Min      I 

an  »  •  •  .  »a   i=i 

1'      m 


2 
a,  d),  (x .  ,y .  )  +...+  a  d>  (x.,y.)  -  f. 
lrl  i   J  i mYm   i'-^i' i 

Pi(xfy) 


for  given  (x,y) .   McLain  gave  results  of  a  number  of  tests 

for  several  choices  of  basis  functions,  d>  .  ,  and  distance 

D 

functions,  p..   The  dp.  consisted  of  the  low  order  monomials, 
1  3  n 

k  t  2  ad 

x  y  ,  and  p.  was  taken  as  d . ,  d .  ,  or  d . e   i  .   The  higher 

■*  i  ill  3 

power  of  d.  and  the  exponential  are  motivated  by  the  desire 
to  place  less  importance  on  distant  data  than  would  be 
accomplished  by  d.  alone. 

The  function  (2.1)  is  not  defined  at  any  data  point, 

(x.,y.)f  by  the  above  definition.   Because  the  weight  for 

2 
the  best  least  squares  approximation  1/p .  (x,y)  — ►  °°  as  (x,y) — > 

(xi,yi)  it  is  clear  that  P[f](x,y) — ►  f ^    as  (x,y)  — >■    (x.  ,yi)  , 
else  the  sum  of  the  squared  errors  would  be  unbounded.   Thus 
we  obtain  a  continuous  approximation  if  we  define  P[f] (x.,y.)  = 
f.,  i  =  1,...,N.   McLain  asserts  that  the  interpolants  are 
infinitely  dif f erentiable,  and  upon  the  assumption  of  non- 
singularity  of  the  coefficient  matrix  for  the  normal  equations 


obtained  from  (2.2),  this  readily  follows. 

McLain  singled  out  the  case  of  the  bivariate  quadratic, 

P[f]  (x,y)  =  a,  +  iLx  +  a^y  +  a.x   +  a"5xy  +  agy   ,  with 

ctd 
weight  function  p.  =  d.e   1   as  working  well  for  a  variety 

of  problems  in  the  sense  that  small  deviations  from  certain 
test  functions  were  observed.   Additional  experimentation  [3] 
has  confirmed  this,  but  even  so  this  method  has  two  serious 
drawbacks:   (1)  The  computational  effort  required  is  large 
since  each  evaluation  requires  the  solution  of  a  least  squares 
problem,  and  (2)  the  method  is  global,  that  is,  the  interpolant 
depends  on  all  data  points  regardless  of  how  far  away  these 
points  are  from  the  point  at  which  the  interpolant  is  being 
evaluated. 

Both  of  the  methods  we  propose  can  be  viewed  as  modifica- 
tions of  the  above  inverse  distance  weighted  least  squares 
interpolant  in  that  they  consist  of  replacing  a.(x,y)  with 
an  approximation  A[a.] (x,y)  which  is  computationally  more 
tractable  than  a.(x,y)  itself.   We  note  that  as  long  as 


A[a.](x.,y.)  =  a.(x.,y.)» 
]    l  2i      j      l  -*i 


i  =  1, . . . ,N, 


the  operator 


m 
Q[f](x,y)  =    I        A[a.]  (x,y)cf)  (x,y) 
j=l      :        ^ 


will  maintain  the  interpolatory  properties  of  P[f].   The  type 


of  approximation  we  use  for  a.,  j  =  l,...,m  is  of  the  form 

_  N  _ 

A[a.](x,y)  =  I      a. (x. ,y.)  W.  (x,y) 
J  i=i   J   1  x         2- 


where 


(2.3)      W±{x.,y.)  =  6ij/  i,  j  =  1,...,N. 


Therefore,  we  may  rewrite  the  above  expression  for  Q  as 


N 
(2.4)      Q[f](x,y)  =  I      Wi(x,y)Qi(x/y) 


i=l 

m 
I 


where  Q.(x,y)  =   £   a  .  (x . , y^) $ •  (x,y)  is  the  inverse  distance 


weighted  least  squares  fit  at  the  point  (x.,y.). 

We  refer  to  the  functions  Q.,  i  =  1,...,N  as  the  nodal 
functions  since  they  are  associated  with  the  nodes  (x.fy.), 
i  =  1,...,N,  respectively.   We  note  that  Q.(x.,y.)  =  f., 
i  =  1,...,N,  and  that  Q.(x,y)  is  a  local  approximation  (near 
(x.,y.))  to  f(x,y),  and  as  such  may  be  expected  to  mimic 
the  shape  of  f(x7y)  provided  distant  points  do  not  influence 
Q. (x,y)  too  much. 

In  order  for  the  interpolant  Q[f]  to  maintain  the  local 
shape  characteristics  of  the  nodal  functions  we  will  require 
certain  properties  for  the  W.  in  addition  to  (2.3).  Specif- 
ically, to  preserve 


»aiau.,y.,  =  ^i(x.,y.) 

1  '  ■*  1      dx  1  J  1 


3x 


and 


3Q[f] /      x    8Qi/     x   •    i 


we  require  that 


3W.  3W. 

(2.5)      g3p  (xj/Yj)  =  ^i(Xjfyj)  =  0,  i,  j  =  1,...,N 


For  our  two  methods,  we  also  propose  the  use  of  bivariate 
quadratics  for  Q.(x,y),  i  =  1,...,N.   Thus,  if   f   itself  is 
a  quadratic  function,   the  function  Q.  will  be  identical  to 
f,  i.e.,  Qi(x,y)  =  f(x,y),i  =  1,...,N. 
Therefore 


N  N 

Q[f](x,y)  =  I    W.  (x,y)Q.  (x,y)  =  f(x,y)  I    W.  (x,y) 
i=l  X       1  i=l 


whenever   f   is  a  quadratic  function.   In  order  to  obtain 
quadratic  precision  of  the  modified  interpolant  Q[f],  we  add 
the  additional  requirement  that 

N 
(2.6)      I    W.  (x,y)  =  1. 
i=l  2- 

The  choice  of  the  distance  function  to  be  used  in 
(2.2)  when  calculating  the  nodal  functions  was  made  on  the  basis 


of  extensive  numerical  testing  [4].   Franke  and  Little  [2, 
p.  112]  proposed 


(2.7) 


1  _  (Ri  -  di>  + 
p.  ~     R.d. 


(R.  -  d.) , 
i    i  + 


R.  -  d.  ,   R.  -  d.  >  0 
l     11     l 

0        R.  -  d.  <  0  , 
l     l 


and  for  appropriate  choice  of  the  values  R.  ,  this  works  well. 


We 


recall  that  Q,  is  the  solution  of  the  inverse  distance  weighted 


least  squares  problem  at  (x,y)  =  (x,  ,y,  ) .   Discussion  of  the  Q, 
is  simplified  (as  is  the  actual  computation)  by  assuming 
that 


Qk(x,y)  =  fk  +  ak2(x-xk)  +  aR3 (y-yk)  +  ak4 (x-xR)   + 


ak5(x-xk)(y-yk)  +  ak6 (y-yk) 


We  then  seek  the  solution  of 


Min 


ak2'*** 'ak6 


N 

I 
i=l 


fv  +  auo(x^-xu)  +..-+  a,c(y.-y.)   -  f. 


1      2 


k2 v  i  ^k 


k6wi  ■'k' 


pi(xk'V 


If 


is  given  by  (2.7),  then  whenever  d.(x,  ,y,)  >  R.  , 

1   K.   K,        1 


the  point  (x.,y.,f.)    has  no  influence  since  the  corresponding 
term  in  the  sum  is  zero.   Thus  Q,  depends  only  on  "nearby" 
points  and  is  therefore  a  local  approximation  to   f.   If  we  now 
use  weight  functions  W,  which  are  nonzero  only  in  some  neigh- 
borhood  of  (x,  ,y,  )  we  will  obtain  a  local  interpolant. 

The  proper  choice  of  the  "radii  of  influence",  R., 
is  critical  to  the  success  of  the  method,  and  we  will  discuss 


this  in  the  context  of  our  first  method  in  the  next  section, 
as  well  as  in  Section  5. 

3.   Two  Methods  for  Interpolation  To  Scattered  Data 

Specifying  the  weight  functions  W.(x,y),  i  =  1,...,N  to 
be  used  in  (2.4)  will  define  an  interpolant.   We  have  found 
that  both  of  the  choices  we  propose  have  very  comparable 
fitting  capabilities,  but  we  feel  that  there  are  situations 
in  which  one  or  the  other  may  be  preferable,  and  so  both  are 
presented. 

3.1   Method  I. 

This  scheme  is  based  upon  a  special  case  of  the  inverse 
distance  weighted  least  squares  interpolant  given  by  (2.1)  - 
(2.2).  If  m  =  1  and  (J),  (x)  =  1,  then  (2.2)  can  be  explicitly 
solved  to  yield 


N   f . 

I      — 
i=l  p  .  (x,y) 


a1(x/Y) 


N   l 

L         2 
i=l  p±   (x,y) 


and  so 


(3.1)     P[f]   _ 


N   f . 

N    i 

1=1  Pi 


This  method  was  first  proposed  by  Shepard  [11] .   Without 
modification  it  does  not  have  very  good  fitting  properties. 
Gordon  and  Wixom  [5]  have  analyzed  this  method  and  have  proposed 
some  interesting  modifications.   They  also  discuss  some  application 
areas  that  are  well  suited  for  this  type  of  interpolant.   The 
fact  that  Shepard 's  method  is  a  special  case  of  inverse  distance 
weighted  least  squares  interpolants  has  been  pointed  out  by 
several  authors,  e.g.  [10].   It  can  easily  be  verified  that  as 
long  as  p.(x.,y.)  =  0,  the  functions 

(3.2)      W.  =     pi 
l 


N   l 

k-i  4 


will  satisfy  the  conditions  of  (2.3),   (2.5),  and  (2.6).   Th 


is 


last  statement  assumes  that  the  distance  functions  p.,  i  =  1,...,N 

are  sufficiently  smooth  so  that  the  derivatives  (2.5)  exist. 

This  is  certainly  the  case  for  our  choice  of  p.,  i  =  1,...,N 

given  by  (2 . 7) . 

We  complete  the  description  of  our  first  method  with  a 

discussion  of  the  selection  of  the  R.  involved  in  the 

definition  of  p . .   While  the  use  of  variable  radii  (i.e.,  R.  f 

R.)  adds  to  the  flexibility  of  the  methods,  we  have  found  that 

for  a  general  purpose  interpolant,  the  selection  of  these  values 

can  quite  often  be  a  nuisance  which  can  be  avoided  by  the  use  of 

a  uniform  value  (R.  =  R  for  all  i) . 

l 

At  this  point  we  emphasize  that  the  distance  functions 
p.  enters  in  two  places:   (1)  definition  of  the  nodal 


10 


functions,  Q,  ,  and  (2)  definition  of  the  weight  functions,  W,. 

It  is  not  necessary  to  use  the  same  radius  of  influence  for 

both  instances,  and  experience  has  shown  it  is  desirable  to  use 

different  values,  say  R  =  R   in  defining  the  nodal  functions, 

and  R  =  R   in  defining  the  weight  functions.   Since  R  denotes 
w  ^         ^  q 

the  radius  of  influence  of  the  data  points  on  the  nodal  functions, 

while  R  denotes  the  radius  of  influence  of  the  nodal  functions 
w 

on  the  interpolant  Q[f] (x,y) ,  it  is  clear  we  should  take  R   < 
R  .   In  order  to  aid  the  naive  user  in  making  reasonable  choices 

q 

of  R  and  R  we  have  found  it  convenient  to  specify  values  of 
q      w  c     j 

N   and  N   and  to  compute  the  radii  of  the  influence  regions 
according  to  the  relationships: 


Rq  ~  2 


V 


(3.3) 


Rw  =  iW 


V  N 


where  D  =  max  d. (x . ,y . ) . 
i/D      J   J 

These  choices  of  R   and  R   eliminate  the  effects  of  scaling 

q      w  ^ 

the  data.   The  values  of  N   and  N   can  be  thought  of  as 

q       w  3 

representing  the  number  of  data  points  anticipated  to  lie  in 

circles  of  radii  R   and  R  ,  respectively.   For  somewhat  uniformly 

q       w '  r  j  j 

distributed  data,  we  have  found  that  a  value  of  N   =18  works 

q 

quite  well.   For  data  that  have  some  regions  which  are 
relatively  sparsely  populated  and  other  regions  where  the 
data  are  comparatively  dense,  or  for  small  sets  (N  <  25) ,  it 
may  be  necessary  to  increase  the  values,  since  the  interpolant 


11 


is  defined  only  on  the  union  of  disks  of  radius  R   centered 


at  the  data 


points  (xi,yi)/  i  =  1,...,N 


We  have  also  found 
2  is  useful.   To 
avoid  problems  when  fewer  than  five  additional  data  points 


that  the  use  of  the  relationship  N  /N 

CJ     Vv 


fall  within  a  distance  R  of  some  given  (x.,y.)/  we  have 

q  ^        i'Ji" 

incorporated   an  automatic  fallback  to  linear  nodal  functions 
in  this  case.   In  general,  we  use  the  singular  value  decomposi- 
tion to  compute  the  coefficients  of  the  nodal  functions,  which 
avoids  possible  nonuniqueness  problems.   It  is  not  usual  for 
either  situation  to  occur,  however. 

We  now  summarize  the  description  of  our  first  method. 

i)    Select  N   and  N   in  order  to  define 
q      w 


(R  -d.) . 

_3 i_t 

R  d. 
q  i 


q    2  W  N 


(R  -d 
w 


R  d 

w 


^  <  •  =  i J¥ 

w     2   W  N 


Default  values  of  N   and  N   are  18  and  9,  respectively 

q         W  '  tr  J 

ii)   For  k  =  1,...,N  solve  the  least  squares  problem: 


min 


a,  .  ,  j  —  2 ,  .  .  .  ,6 


N 


/Li     P?(*v' 


i=l 
i^k 


i^'V 


fk+ak2(xi-xk)+ak3(yi-yk) 


-i  2 


+ak4  (xi-xkJ  2+aks  (xi-xk}  ^i-ykJ  +ak6  ^r^]  2-f  i 
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to  yield  a,  .  ,  j  =  2,..., 6. 
iii)  Define 


Qk(x,y)  =  fk  4-  ak2(x-xk)  +  a^  (y-yk) 

+  ak4(x-xk)2  +  ak5(x-xk)  (y-yk)  4-  aR6  (y^)  2 
and  compute 

N  Qk(x,y) 

kLp^(X/y) 
D[f] (x,y)  = 


N      x 
k=l  Pk(x^} 


3.2   Method  II. 

Our  second  method  makes  use  of  a  triangulation  of  the 
data  V.  =  (x.,y.)  i  =  1,...,N  in  order  to  define  the  weight 
functions  W.,  i  =  1,...,N.   We  use  an  algorithm  which  triangulates 
the  convex  hull  based  upon  the  min-max  angle  criterion  as 
described  by  Lawson  [6].   A  FORTRAN  program  which  implements 
this  algorithm  is  available  as  part  of  [1] .   Alternatively,  if 
a  triangulation  is  in  existence  for  other  purposes,  it  can  be 
used. 

Each  W.  will  be  a  globally  defined  C   function  with  support 

S.  =   U       T-i,i'  where  T.,,  denotes  the  triangle  with 
1    jkl  £  M.   DK1         DK1 

vertices  V.,  V,  ,  and  V,  and  M.  =  {jkl:  T.,  ,  is  a  triangle 
j         k'       1       l     J      jkl  ^ 

with  vertex  V . } . 
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Figure  1:   S 


We  first  define  W.  and  its  first  order  partial  derivatives  on 

E  =    U    e,  .  where  e,  .  represents  the  edge  with  vertices 

kjeN    K:         K:> 
J      e 


V,  and  V.  and  N   =  {kj:  V,  V .   is  an  edge  of  the  triangulation} 

Following  this,  we  incorporate  a  blending  method  for  triangles 

to  extend  the  definition  to  the  interior  of  each  triangle  of 

the  triangulation.   Let  e. .  be  an  edge  contained  in  S.  with 

lj         ^  l 

V.  as  an  endpoint.   As  a  univariate  function  along  this  edge, 
W.  must  satisfy  four  conditions  imposed  by  (2.3)  and  (2.5). 
Namely 


and 


W. (V.)  =  1,  W. (V.)  =  0, 
li         i   j 


aw.  9w. 

i      i  8x    i'    vjrn  n  9v    i 


aw.  aw. 
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These  conditions  can  be  satisfied  by  a  cubic  polynomial  and 
so  we  define 

Wi((l-t)Vi  +  tVj)   =    (1-t)2  (2t+l),  0  <  t  <  1. 

On  all  other  edges  we  define  W.  to  be  zero.   In  order  to 
maintain  continuity  of  the  first  order  derivatives  across 
edges /  it  is  convenient  to  specify  the  derivatives  normal  to 
an  edge  as  a  linear  function  along  the  edge.   That  is 


3W.  3W. 

(yj-ri>33r((1-t)vi +  tvj»  -  (xrxi»F((1't,vi +  tvj' 

aw.  aw. 

=  (l-tUfy.-y.^V,)  -  (Xj-xi)wi(V.)] 


9W.  9W. 

+  t[(y.-y.)-5-^(V.)  -  (x.-x.)tt— ^(V.)  ] 
Jj    ;i3x    j '      j      i  3y    j 


In  light  of  (2.5),  this  means  that  the  normal  derivatives  will 
all  be  identically  zero.   This  completes  the  description  of 
the  edge  information  for  W.. 

In  order  to  extend  the  definition  of  W.  to  the  interior 

i 

of  each  triangle,  we  use  an  interpolation  method  [9]  which 

will  assume  predescribed  position  and  slope  on  the  entire 

boundary  of  a  triangle  domain.   After  substituting  the  edge 

information  into  this  triangular  blending  method,  we  find 

that  for  (x,y)  e  T. ..  c  s.,  the  weight  functions  have  the 
1  ljk     i         ^ 


form 


W.  (x,y)  =  b2(3-2b. ) 
l   J  l      l 


(3.4) 


+  3 


b2b.b, 
1  J  k 


b.b  .+b.b1  +b  .b. 
l  j      i  k  j    k 


e . 
l 


+   e 


lejll2 


where  b.,  b.,  b,  are  the  barycentric  coordinates  of  (x,y] 

1  K 


with  respect  to  the  triangle  T. .,  and 

c  ^  ljk 


n 


n  =  i ,  j  or  k 


represents  the  length  of  the  edge  opposite  V  ,  n  =  i,  j  or  k. 
The  barycentric  (area)  coordinates  are  given  by  the  equations 
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ejll    2    +     llekll  2    -   ll^ll  2 


-      e, 


2~| 


x   =   b.x.    +   b.x.    +   b,  x, 
li  j    j  k   k 


(3.5) 


k*k 


y   =   b.y.    +   b.y.    +   b,y 

i    !  ]    ] 


1   =   b .    +   b .    +   b,     . 
13k 


We  can  now  note  that  on  an  arbitrary  triangle  T. .,  the  only 
weights  which  are  nonzero  are 
final  interpolant  is  given  by 


weights  which  are  nonzero  are  W.,  W.  and  W,  .   Therefore,  the 
*  ljk 


(3.6)      G[f](x,y)  =  Wi(x/y)Qi(x,y)  +  W^ (x,y) Q^ (x ,y) 

+  Wk(x,y)QR(x,y)  ,  (x,y)  £  Tijk- 

We  also  note  that  W.  +  W.  +  W.  =  1,  so  that  (2.6)  is  satisfied. 

ljk 

We  now  summarize  this  method. 

1)    Define  the  nodal  functions  Q.  ,  k  =  1,...,N 

as  in  Method  I, 

ii)    Form  a  triangulation  of  the  points  V.  =  (x. ,y.) , 

i  =  1, . . . ,N, 

iii)    Given  (x,y)  determine  the  vertices  V.,  V.  and  V, 

1    J       K 

of  the  triangle  that  contains  this  point  and 

compute  G[f](xfy)  according  to  (3.6)  using  (3.4)  and  (3.5) 
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We  note  that  this  method  is  very  similar  to  that  proposed 
by  McLain  [8].   Our  weight  functions  W.,  i  =  1,...,N  arise  in 
a  natural  manner,  however,  and  use  of  the  distance  function 
given  in  (2.7)  is  an  improvement  since  it  is  generally  the  case 
that  shape  characteristics  of  the  nodal  functions  in  relation 
to   f   are  adversely  influenced  when   global  approximations 
are  used. 

4 .    Examples 

In  order  to  illustrate  the  performance  of  the  two  methods, 
we  include  some  examples.   These  examples  are  only  a  few  of 
many  upon  which  our  conclusions  are  based,  but  are  representative 
of  the  methods'  approximation  properties.   The  first  group  of 
examples  utilizes  ordinates  from  the  function 


f(x,y)  =  .75EXP  [-  (9x-2)2+(9y-2)2j  +  <75EXp 


-  .2EXP    j"-(9x-4)2  -  (9y-7)2J+  . 


5EXP 


(9x+l)2    (9y+l) 
49     ~    10 


(9x-7)2+(9y-3)2 
4 


A  perspective  plot  of  this  surface,  viewed  from  the  first 
quadrant  at  an  angle  of  30°  from  the  x-axis  is  given  in  figures 
3a  and  4a.   The  function  was  chosen  so  as  to  present  a  variety 
of  behavior  in  a  single  surface.   The  maximum  value  of  the 
function  is  approximately  1.22  near  the  point  (.22,  .22),  while 
the  minimum  value  is  approximately  .004  near  the  point  (.47,  .78). 
Three  sets  of  data  (x.,y.),  i  =  1,...,N  were  used. 
Set  1:   This  set  consists  of  100  points  generated  by  a 
pseudorandom  number  generator,  one  point  in  each  subsquare 
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(c)    Data    Set    3 


Figure    2. 


18 


of  side  1/9  centered  at  (i/9,j/9),  i,  j  =  0,1,. ..,9.   These 
points  are  shown  in  Figure  2a. 

Set  2:   This  set  of  33  points  was  selected  manually  to  have 
regions  of  varying  density.   These  points  are  shown  in 
Figure  2b. 

Set  3:   This  set  of  25  points  was  selected  manually  to  yield 
a  somewhat  uniform  coverage  of  the  unit  square,  and  are  similar 
to  a  set  appearing  in  [8] .   These  points  are  shown  in  Figure  2c 

The  interpolants  were  evaluated  on  a  uniform  grid  of 
33  x  33  points  in  the  unit  square.   The  resulting  surfaces 
are  shown  in  Figures  3  and  4.   Table  1  contains  the  maximum, 
mean,  and  RMS  deviations  at  these  points.   For  Method  II  a 
few  of  the  display  points  lie  outside  the  convex  hull.   For 
plotting  purposes,  these  were  set  to  zero  and  were  omitted 
in  the  calculations  leading  to  the  errors  of  Table  1. 


Case 

Max 

Mean 

RMS 

1.1 

.  0573 

.0079 

.0128 

1.2 

.1844 

.0340 

.0478 

1.3 

.1584 

.0353 

.0486 

II. 1 

.0481 

.0072 

.0113 

II. 2 

.1501 

.0326 

.0455 

II. 3 

.1535 

.0349 

.0475 

Table  1.   Errors 


19 


Figures  3b  and  4b  show  that  both  methods  reproduce  the 
surface  quite  well  for  data  set  1.   The  main  defect  appears 
to  be  near  the  higher  peak,  particularly  for  Method  I.   As 
can  be  seen  from  the  disposition  of  the  points  in  Figure  2a, 
this  is  in  a  region  where  a  relatively  large  gap  between  points 
occurs.   A  similar  occurrence  accounts  for  the  depressed  area 
behind  the  dip  in  both  figures  3b  and  4b. 

Figures  3c  and  4c  are  quite  similar  and  both  have  noticeable 
defects  when  compared  to  the  test  surface.   In  particular,  the 
dip  is  completely  missed  because  the  data  points  fail  to  define 
it.   The  appearance  of  a  crease  in  Figure  4c  near  the  right 
rear  edge  is  due  to  the  occurrence  of  a  long  thin  triangle 
along  that  edge  which  causes  blending  of  two  nodal  functions 
from  points  relatively  far  apart  in  the  definition  of  the 
interpolant  for  that  triangle.   When  these  two  nodal  functions 
differ  significantly,  as  they  do  here,  rapid  variations  can 
occur  across  the  narrow  part  of  the  triangle. 

Figures  3d  and  4d  appear  to  be  less  alike  than  they  actually 
are  because  values  outside  the  convex  hull  have  been  set  to 
zero  in  Figure  4d.   The  most  significant  difference  between 
the  two  is  near  the  left  rear  edge,  where  Figure  3d  shows 
the  surface  (apparently)  dipping  down  rather  rapidly,  while 
Figure  4d  shows  a  near  crease  similar  to  that  in  Figure  4c. 
Because  of  a  data  point  near  (.47,  .78),  the  dip  is  partially 
defined  in  this  case,  but  since  there  are  no  other  nearby 
points  it  is  extended  over  a  much  wider  area  than  in  the 
test  surface. 
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(a)  Test  surface 


(b)  Interpolant  for  Data  Set  1 


(c)  Interpolant  for  Data  Set  2 


(d)  Interpolant  for  Data  Set  3 


Figure  3:  Method  I 
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(a)  Test  surface 


(b)  Interpolant  for  Data  Set  1 


c)  Interpolant  for  Data  Set  2 


(d)  Interpolant  for  Data  Set  3 


Figure  4 :  Method  II 


Another  set  of  data  obtained  from  Akima  [1] ,  which  arises 

in  a  study  of  waveform  distortion,  is  given  in  Table  2.   We 

use  this  set  to  illustrate  the  effects  of  varying  the  parameters 

N   and  N  .   The  results  are  shown  in  Figure  5. 
q      w  3 

Figures  5a  and  5c  appear  to  be  difficult  to  choose  between. 
Very  slight  differences  can  be  observed  along  the  front  edge. 
Figure  5e  is  definitely  less  desirable  than  5a  or  5c  because  of 
more  undulations  near  the  front  edge  and  a  higher  peak  at  the 
right  rear.   Extensive  testing  has  shown  that  Method  I  is 
fairly  stable  for  values  of  N   and  N   in  the  ranges  given  here. 

Figure  5d  appears  to  be  the  more  desirable  surface  among 
5b,  5d,  and  5f.   Figure  5b  shows  a  small  defect  near  the  right 
front  edge,  while  the  choice  between  5d  and  5f  is  less  obvious. 
All  three  show  the  characteristic  defect  over  the  long  slim 
triangle  at  the  right  rear  edge,  allowing  the  surface  to  dip 
out  of  sight  there.   It  should  be  pointed  out  that  there   is 
no  known  "parent"  surface  here,  and  in  fact  that  behavior  may 
be  correct,  although  that  part  of  the  interpolant  is  not 
pleasing  because  of  the  very  rapid  changes  which  occur. 

5.    Conclusions  and  Recommendations 

The  two  schemes  discussed  here  have  been  found  to  be 
capable  of  generating  reasonable  interpolation  functions  in 
a  variety  of  cases.   A  number  of  comments  are  appropriate, 
however. 

First,  we  have  restricted  ourselves  here  to  the 
discussion  of  local  interpolants.   Local  methods  are  necessary 
for  very  large  sets  of  data,  but  in  exploring  their  properties 
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i 

x  . 

l 

*i 

z  . 

i 

i 

x  . 

l 

yi 

z  . 

l 

1 

11.16 

1.24 

22.15 

26 

3.22 

16.78 

39.93 

2 

24.20 

16.23 

2.83 

27 

0.00 

0.00 

58.20 

3 

12.85 

3.06 

22.11 

28 

9.66 

20.00 

4.73 

4 

19.85 

10.72 

7.97 

29 

2.56 

3.02 

50.55 

5 

10.  35 

4.11 

22.33 

30 

5.22 

14.66 

40.  36 

6 

24.67 

2.40 

10.25 

31 

11.77 

10.47 

13.62 

7 

19.72 

1.39 

16.83 

32 

17.25 

19.57 

6.43 

8 

15.91 

7.74 

15.30 

33 

15.10 

17.19 

12.  57 

9 

0.00 

20.00 

34.60 

34 

25.00 

3.87 

8.  74 

10 

20.  87 

20.00 

5.74 

35 

12.13 

10.79 

13.71 

11 

6.71 

6.27 

30.97 

36 

25.00 

0.00 

12.00 

12 

3.45 

12.78 

41.24 

37 

22.33 

6.21 

10.25 

13 

19.99 

4.62 

14.72 

38 

11.52 

8.  53 

15.74 

14 

14.26 

17.  87 

10.74 

39 

14.59 

8.71 

14.81 

15 

10.28 

15.16 

21.59 

40 

15.20 

0.00 

21.60 

16 

4.51 

20.00 

15.61 

41 

7.54 

10.69 

19.31 

17 

17.43 

3.46 

18.60 

42 

5.23 

10.72 

26.  50 

18 

22.80 

12.39 

5.47 

43 

17.32 

13.78 

12.11 

19 

0.00 

4.48 

61.77 

44 

2.14 

15.03 

53.10 

20 

7.58 

1.98 

29.87 

45  - 

0.51 

8.37 

49.43 

21 

16.70 

19^65 

6.31 

46 

22.69 

19.63 

3.25 

22 

6.08 

4.58 

35.74 

47 

25.00 

20.00 

0.60 

23 

1.99 

5.60 

51.81 

48 

5.47 

17.13 

28.  63 

24 

25.00 

11.87 

4.40 

49 

21.67 

14.36 

5.  52 

25 

14.90 

3.12 

21.70 

50 

3.31 

0.13 

44.08 

Table    2. 
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(a)  Method  I:  N  =12  ,N  =6 
q      w 


(b)  Method  II:  N  =12 

q 


(c)  Method  I:  N  =18,  N  =9 
q      w 


(d)  Method  II:  N  =18 

q 


Figure  5. 
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(e)  Method  I:  N  =  24 ,N  =12 
q      w 


(f)  Method  II:  N  =24 

Si 


Figure  5.  (continued) 
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one  need  not  consider  large  sets  because  they  are  local. 
We  have  only  considered  interpolation  methods  here,  although 
we  recognize  that  often  it  may  be  more  appropriate  to  use 
approximation  methods  which  smooth  the  data  in  some  sense. 
While  we  have  not  investigated  this  possibility,  it  is  clear 
that  replacement  of  the  nodal  functions  Q,  (x,y)  with  local 
smoothing  functions  rather  than  ones  which  take  on  the  value 
f,  will  lead  to  a  smoothing  approximation  which  is  local. 

The  use  of  the  same  radius  of  influence  for  each  point, 
and  the  calculation  of  that  radius  from  the  diameter  of  the 
point  set  is  strictly  a  matter  of  convenience  for  the  user. 
One  could  argue  that  this  device  makes  the  methods  global 
since  addition  of  a  point  will  change  the  radius  of  influence 
and  hence  the  entire  interpolant.   For  this  reason  we  made  the 
computation  of  the  radius  of  influence  an  option  (although 
it  is  the  default  option) .   The  use  of  different  radii  of 
influence  could  be  a  necessary  and  desirable  feature  when  the 
density  of  points  varies  drastically  over  the  point  set.   Our 
experience  indicates  that  one  should  probably  choose  radii  so 
that  the  disks  associated  with  a  point  contain  approximately  18 
points  for  defining  the  nodal  functions  (Methods  I  and  II)  and 
about  9  points  for  defining  the  weight  functions  (Method  I) . 

Regarding  the  choice  between  Method  I  and  Method  II, 
we  make  the  following  comments:   Method  I  has  the  advantage  of 
simplicity.   While  Method  II  requires  a  triangulation  (and  the 
machinery  for  generating  it  if  it  is  not  already  available) , 
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and  thus  more  auxiliary  storage,  it  is  considerably  faster 
since  each  evaluation  involves  only  three  nodal  functions, 


while  Method  I  typically  involves  about  9  (N   =9)  nodal 


functions.   Method  II  also  has  the  disadvantages  noted  in 
the  examples  when  long  thin  triangles  occur.   Method  II  is 
not  readily  extended  to  more  than  two  independent  variables, 
as  is  Method  I.   Nonetheless,  for  certain  applications,  Method 
II  may  be  the  appropriate  choice,  particularly  if  a  triangula- 
tion  is  already  in  existence. 

In  conclusion,  we  note  that  all  local  methods  involve 
some  ad  hoc  assumptions  and/or  parameters  to  be  specified 
by  the  user.   The  methods  we  propose  have  endured  through 
extensive  testing  of  their  fitting  properties  and  appropriate 
values  for  their  parameters.   We  feel  they  will  perform  quite 
adequately  in  a  variety  of  situations.   Nonetheless,  we 
recognize  that  selection  of  a  suitable  interpolant  is  a  sub- 
jective matter  even  in  the  case  of  one  independent  variable, 
and  thus  the  choice  of  an  interpolant  ultimately  rests  with  the 
user. 

FORTRAN  programs  which  implement  Methods  I  and  II  are 
available  on  request  from  the  authors. 
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of  this  paper. 
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